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The performance of single-phase storage and expulsion sys- 
tems is strongly affected by temperature variations within 
the stored cryogen which are generated during heat trans- 
fer. Peculiar operating responses are indicated by spon- 
taneous changes in fluid pressure which accompany "g" level 
changes, increased heater surface temperature, and dura- 
tions of pressure cycles which differ considerably from 
that which is computed for an isothermal cryogen. The non- 
isothermal characteristics are predicted with a numerical 
model which includes the simultaneous solution of the time 
dependent conservation equations of mass, energy, and mo- 
mentum in two space dimensions of Cartesian coordinates for 
boundary conditions which approximate those of the flight 
cryogenic system. The methodology of the numerical method 
and some comparisons between the predictions and the Apollo 
12 flight data are included. 
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PRESSURIZED EXPULSION OF NONISOTHERMAL 

SINGLE-PHASE CRYOGEN 

Clifford K. Forester 
Research Engineer 
The Boeing Company 
Seattle, Washington 


The single-phase storage concept has been employed in the 
design of the Gemini and Apollo atmosphere and fuel cell 
supply equipment. In the course of using these systems, 
some peculiar performance characteristics have been ob- 
served during flight operations. For example, shortly 
after the launch of Gemini II and Apollo 12, abrupt de- 
creases in the oxygen storage pressure in excess of 100 psi 
have been noted. Additionally, the observed time required 
to complete a pressurization cycle has been both substan- 
tially shorter and longer than the duration predicted from 
calculations which assume the stored cryogen is isothermal. 

In this paper, the methodology of a numerical model is pre- 
sented which, for the comparisons this far made with Apollo 
12 and 14 flight data, appears to predict the observed pe- 
culiarities. Several considerations comprise the numerical 
modelx 

1) the methodology of the General Elliptic Method - GEM 

2) the 1/2 box model in which the GEM is applied 

3) the method of accounting for the cryogenic container 
elasticity 

4) the method of treating the gas trapped in the external 
volume 

The first two of these are discussed in order after the 
introduction. The mathematical derivation of the last two 
are provided in Appendix A. These considerations are fol- 
lowed by a discussion of the physical characteristics of 
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the Apollo oxygen system and a brief discussion of the 
flight data reduction considerations. Finally, a discus- 
sion is provided in which the comparison between Apollo 12 
flight data and the 1/2 box model predictions for two pro- 
blems is presented. 

The Methodology of the General Elliptic Method (GEM) 

Because of the length, no attempt is provided here to docu- 
ment the literature now available on numerical algorithms 
which model the conservation equations of mass, energy, and 
momentum in time and space. A review of this literature 
at the end of 1967 showed that the methods generally ap- 
plied either to incompressible flow or high speed compres- 
sible flow. These algorithms are restrictive enough so 
that they eliminate from practical analysis one of the 
most intriging aspects of the problems at hand which is to 
compute the cryogen pressure history. As a result, a new 
algorithm, the General Elliptic Method (GEM) , was developed 
by the author as a part-time effort and was documented in 
reference (1) . Reference (2) includes a literature survey 
of those algorithms which depend upon an elliptical equa- 
tion and some refinements of the GEM which was presented 
in reference (1) . However, the GEM had to be modified to 
suit the special requirements which the problems at hand 
establish. Since these modifications are an intricate 
part of the computational method, a complete description 
of the methodology is necessary and is illustrated subse- 
quently. 



The conservation equations of mass, energy, and momentum 
for a Newtonian fluid may be written in the x- and y- 
directions of Cartesian coordinates respectively as 


ap . 2>Ci°a) . a(pvQ _ 0 
at ax 
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at + ax 

6(<qu) a(pa 2 ) 
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where 



and where p is the mass density, u is the x-direction 
velocity, v is the y- direction velocity, e is the speci- 
fic internal energy, h is the enthalpy, k is the thermal 
conductivity, T is the temperature, P is the pressure, g x 
is the x- direction acceleration component, gy is the y- 
direction acceleration component, and ytt is the kinematic 
viscosity. 


Since 

e « c (P, P) and h ■ c + P/p 


where 

v p ac 
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■ ■ -MW, 


and the product of $ and 0 
squared of the fluid. 


( 4> o ) 
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is the speed of sound 
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Because the Mach number is less than about 10“ 3 for the 
problems to be considered, only pressure waves of the 
acoustic variety will dominate during the various flow pro- 
cesses of interest. But since, for the purposes here, 
modeling acoustic waves is unnecessary, and since each 
acoustic wave only subtly modifies the density field as it 
passes by any point of interest, the local time derivative 
of pressure (bP/dt) may be replaced with the time deriva- 
tive of the global pressure (P g ) or ( a P g / & t) -► p . ( Note 

however that this simplification does not eliminate spatial 
gradients of pressure from the momentum equations.) Thus 
equation (2) may be simplified to read 




(3) 


(4) 


Now the mass rates which pass normal to and through the 
center of the sides of a cell are defined as 

X s and y ■ <°\rAy 

where A x and Ay are the areas of the cell faces perpendi- 
cular to the x and y axes respectively. Equations (1, 3 
and 4) can be combined as 

A'° ,i ax ,i aY. _ n 

at lt x 2>X Au au ° (5a) 




at 




£ + A* P » A a o- a 


at 2>\j i* ~*T N =A 






(5c) 


( 5d) 


Now (5a) - (5d) are transformed into a discrete set of 
simultaneous algebraic equations with the use of a discrete 
grid network. 


It should be noted that two types of grid networks are in 
general use: 
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1) a mathematical grid in which all dependent variables 
are defined at common points and 

2) a physical grid in which the state and transport 
properties are defined in the centers of control 
volumes and the velocities are defined normal to and 
in the center of the faces which bound the volumes. 

The physical grid possesses three important advantages: 

1) a unique definition of local mass conservation is 
specified. 

2) The resolution of thermal and velocity boundary 
layers are improved for the same grid density since 
the temperature and velocity points are Ax/2 or Ay/2 
from solid boundaries rather than the usual A x or 
A y for the mathematical grid. This is achieved 
without the usual complication that nonlinear grid 
spacing s involve. 

3) The procedures for treating mathematical singulari- 
ties (e.g., the center of a polar coordinate system) 
are not ambiguous. 

For these reasons the physical grid is employed in the fol- 
lowing discrete formulation. Figure (1) shows a typical 
computational cell which is imbedded in an array of such 
cells which altogether comprise the volume of the entire 
region of interest. For identification purposes, the sub- 
scripts i and j are used with the cell center variables 
P, p , T, h, e, k, and p. to denote their relative location 
in the complete cellular array. The whole integers, i and 
j, are counted with the increasing x- and y- directions, 
respectively. The velocities and mass rates which are de- 
fined at the cell sides are denoted by half integers and 
are counted similarly to the cell center values. Both the 
integer and half integer nodal points are defined to be 
separated by a cell width { A x) and a cell height (Ay). 
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The difference approximations of equations (5a) - (5d) for 
a uniform two dimensional rectangular grid with an arbitra 
ry depth in the third dimension, are 





= p. . + D 


t+l At 
lj V 


(6a) 


<- c vj + *i) °M 


\ At 

‘ V 


energy equation 



(6b) 

(6c) 

(6d) 

(6e) 


where 


oVj' + 

and V = A x Ax= A y Ay(the cell volume) and the superscripts 
t + 1 and t + 6 indicate an evaluation at the time t + A t 
and t + 6 (At) respectively ( 6 = 1/2, 0). For 6 = 1/2, 
i j , 0 i j , Pg and the terms of the convection of heat 
are nearly time centered. Such an alternative has certain 
advantages which are provided together with the methodology 
for implementation in reference (2) . 

D^-; is the net added mass rate to an cell, E^j is the net 
added heat rate to an cell by conduction and convection. 
HXi+ 1 / 2 j and HYij+ 1/2 are the viscous, body force and the 
convection terms of momentum in the x- and y- directions 
respectively. The difference expressions for these terms 
are given in references (1) and (2) and are omitted be- 
cause they are not essential to the illustration of the 
methodology of the GEM. 


58 



( 7 ) 


(6b) may be revised as 
itl • . f vt 

D '-i = % ^ 

Since summed over the total volume is equal to net mass 

rate inflow at the boundaries, and since P g is a constant 
over the total volume, (7) may be integrated over the total 
volume to yield 

p< /net mass rate outflow \ 

' at the b oundaries *_ ^ 






4 

where 23 is the sum over all active cells. 


Now (6a) and (7) may be combined to form 



t At 
+ V 




(4>e) 


ij 



(9) 


Equations (6c, 8, and 9) apply only to a container whose 
walls are rigid. These equations may be modified in two 
ways to accomodate the elasticity effect of the container 
wall and the external volume. The first method involves 
simply the multiplication of (8) by (EF) where 

tF * 

where 


E5F =. 


3r0 ~cr) 
2 bE. 



Pg ft 


a 


( 10 ) 


which is the elastic factor derived in Appendix A. { p<$>& ) c 
is evaluated at the average stored fluid density. ( /OV)p 
is the product of the external plumbing volume and the 
density of the fluid in the volume, n is the inverse of 
the polytropic coefficient and is unity for an isothermal 
process. This method of modifying P g for the elastic 
effect is valid provided the ratio of the volume occupied 
by the thermal boundary layer to the total volume is small. 
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A general method of incorporating the elastic factor in- 
volves modifying (8) as 

_.p»t <net mass rate outflow » 

Z-Aw'ij at the boundaries / 


*" »*>] 
and (9) may be modified as 

P* e / o t . - — hh-. 
r ij * 1 ) V 


( 11 ) 


(12) 


where 


For convenience, equation (6c) is renumbered as 

tvl ^ . • . 

- P# + At tPjg) 


(13) 


Equations (11, 12, and 13) are the discrete conservation 
equations of mass and energy which must be solved simul- 
taneously with the momentum balances (6d) and (6e) . To 
implement this, equations (6d and 6e) are discretely dif- 
ferentiated with respect to x and y respectively and com- 
bine to form 


D ‘j' * A] * - H'^tj > 


(14) 


where 
-ttl * 


Kr *s (p Hr zp u * w ,i+ ^vi 


Mi 


and 


*'A\ * - *»^j - 

The Dij^ term is eliminated from (14) with the aid of (6a) 
and (12) to yield 


0 = + HH*j + At - m\^) 


(15) 
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Equation (15) is an elliptic difference equation which may 
be solved by direct inversion, direct Fourier methods (see 
reference 3) , and iterative methods. For generality and 
simplicity, the Liebmann iterative method is employed 
which is used in the MAC method of Harlow and Welch 
(reference 4) . NIP is defined as the number of iterations 
used to solve equation (15) . The pressures which result 
from the approximate solution of (15) are used in the mo- 
mentum balances to obtain updated x— and y— direction mass 
rates in the active cells. With the use of the thermo- 
dynamic relations, the updated values of density and global 
pressure from equations (12 and 13) are used to find up- 
dated values of temperature, enthalpy, thermal conductivity, 
viscosity, 0 and $ for all active cells. 

The sequence of computations for the General Elliptic 
Method (GEM) is summarized asi 

1) Prescribe the initial values of the mass rates 

(Xi+i/2 j , Yij+i/2) , the global pressure (P g ) , and 
the density (/Oij) for a11 active cells including 
the border values. Compute the velocities (u^g and 
vjLj) for all the field and border locations using the 
mass rate definitions. Evaluate the exterior velo- 
cities for the no-slip and free-slip conditions. Find 
Tij, hij, ( 0 4> ) ij , and 0 i j with P g and p ij 

from the thermodynamic and transport relations. 

These relations may be defined in terms of equations 
of state or tabular data. For the simulations re- 
ported, the thermodynamic properties are tabulated 

at three pressure levels; 850 psia, 900 psia, and 
950 psia. Linear interpolation with density and 
pressure are used to find the other desired thermo- 
dynamic properties. The transport properties are 
determined from a linear interpolation with tempera- 
ture in single arrays of thermal conductivity and 
viscosity at 900 psia. 

2) Evaluate HYij+ 1 / 2 , ®ij , ^ij, (EF) , P g and 

HHij for all active cells from the existing field 
values. 
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3) 


Compute p and Pg t x from (12) and (13). With 
these values find T*^, h t t^/ k /*■ t iT# ( 0 <$> ) 


and 9 from the linear^^ble searche^. 

Compute the border values of X^ + ^/2j and Y ij+l/2 
from the expulsion rate function. 


4) Solve (15) approximately with NIP iterations for 

for all cells inside the border and immediately 
outside the border. The pressure P^j in cells out- 
side and adjoining the border are obtained directly 
from the momentum balances (6d) and (6e) and are 
recomputed each sweep of the field for new tentative 
values. Note that the border values of HXj[+i/2j 
and HYij+1/2 used in (15) , (6d) , and (6e) must be 

the same but may be any desired value. It is best 
to set these to zero. 


5) Use the values of pressure (Pj_j) from step 4 in 

(6d) (6e) to ob^jn the remaining unknown values 

of x i+l/2j and Y ij+l/2 respectively. 

6) The updated values then become the current values at 
time (t) and step 2 through step 5 may be repeated 
until the time (t) has been advanced to some value 
of interest. 


Some discussion is now devoted to the various boundary con- 
ditions which are required to define all derivatives and 
quantities normal to the boundary. All mass rates normal 
to the border must be set to zero except where expulsion 
is specified. The velocities are computed from the mass 
rates, the normal cross sectional area through which the 
mass flows, and the average density of respective adjoining 
cells. Where mass rates are nonzero at the border, the 
velocity calculation requires a prescription of the density 
in the external receiver cell and is set equal to that in 
the supplier cell. This velocity enters into two viscous 
terms. Two other viscous terms require the prescription 
of velocities parallel to the boundary in all external 
cells. Where no-slip boundaries are required, the external 
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velocity is set equal and opposite to the values in the ac- 
tive adjoining cell. Where free-slip boundary velocities 
are required, the external velocity is set equal to the 
values in the active adjoining cell. These velocities in 
external cells are computed after all the velocities have 
been computed which reside inside and on the border. The 
HX i+l/2j and HYij+i /2 terms required by the momentum 
balances on the border may be set to any value but for con- 
venience are set to zero. This follows from the fact that 
since the mass rates at t+1 and t are prescribed on the 
border, the momentum balances are not required to generate 
these values at the border • For these to balance then it 
is only necessary that the values of HX^+l/2j and HY ij+l/2 
defined on the border which are used in the elliptic 
equation (15) be the same as those used in the border mom- 
entum balances to find the unknown external cell pressure 
which is required by equation (15) during the iterative 
solution. 


The heat conduction terms involve gradients of temperature 
normal to the sides of a cell. Cell sides at which the 
heat rate is prescribed may be used to calculate these 
gradients. For generality, it is assumed that any cell 
face at which the heating rate is prescribed has an arbi- 
trary thermal capacitance. The cell face temperature is 
determined with an implicit energy balance on a cell face 
of interest and for a left cell face is 


QX- I . m 




Ax. 

AX 


(16) 


where Qx^ + i/2j i- 3 the prescribed cell face heat rate, 
HCXi=i/2j is the cell face thermal capacitance, lWXi-l/2j 
is the face temperature, k^_3y 2 j the thermal conduc- 
tivity at temperature T wx i_i/2j and T ij the tem- 

perature. Equation (6) may be solved for the unknown 
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t+1 


TWX i-l/2j *>Y algebraic manipulation and yields 


twx 


-V' 




( 17 ) 


Similar expressions to equation (17) may be written for any 
other of the remaining cell faces. It should be noted that 
if (16) were not implicit in the heat conduction term, the 
time step would be restricted to 


At « 


HCX; 


il. 




(18) 


which could be substantially more restrictive than that 
required by GEM. (A discussion of GEM stability conditions 
is given below) . Since the computational efficiency of 
the algorithm is directly related to A t, equation (18) 
could impose severe increases in computing time. This is 
especially important since practical problems frequently 
require hours of computer time. For this reason, equation 
(17) is a preferred form over the explicit alternatives. 
This situation offers a fine example of the reason why 
implicit schemes are preferred or essential compared to 
explicit schemes whether the model of the boundary condi- 
tions or the field conservation equations are concerned. 
Reference (2) has some additional details with regard to 
explicit modeling of the conservation equations. 


Note that for HCXj, . an^X i _ 1/ / 2 j e< 3£ al to zero, equa- 

tion (17) reduces to X D TWXi-l/2j = Tij which is the 
boundary condition for an ideally insulated cell face. 
Where the biased differences for the convective terms, 
reference (2) , are employed, the fundamental stability 
conditions for GEM are: 

(i9) 

where 

A = (k/ p cp) max 

B = 0, 1, 2 for the x- direction, and 



( 20 ) 


C = 0, 
and 
I 

At 


1, 2 for the y-direction 




C; 1 


0 + zfrWx] 


v Ai? + Ay 

B and C depend upon the number of cell faces which have 
prescribed heat rates in equation (19). For " g" spikes, 
CAX, 


At 3 < 




ZjS 


(21) 

are nearly 


is used to insure that equations (19) and (20) 
satisfied during a "g" spike. For time integration accur- 
acy, equation (21) is modified to read 

n a « a » i A 

( 22 ) 


At 4 = 0.03 






In the calculations, equations (19) and (20) multiplied by 
about a factor of 1.5 in order to ensure small truncation 
errors in time while the "g" level is constant. Equations 
(19) and (20) are additionally superceded by the conditions 
If > 0 


then 

At = Min. ( A t]_, A t 2 , 30) 
If At > 1.3 At t-l/2 
then 

A t = 1.3 A t t-1 / 2 


(23a) 


(23b) 


where 

At^ - ^ is the previous time step in seconds. 

Equation (23a) has a thirty second value shown. It is used 
to minimize the truncation errors in time for pressuriza- 
tion strokes which are about four minutes long and which 
are simulated with coarse meshes. 


Equation (23b) is used to limit the rate change of time 
step increases in order to promote time integration accur- 
acy especially after "g* spikes or abrupt changes in the 
heating rate. 
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If there is a "g" spike, then 

A t “ Min ( At, A t 4 ) (23c) 

is imposed for ten time steps. 

There are a number of integral relations which provide 
checks on the accuracy or correctness of the various 
local calculations and these are now discussed. The sum 
over the volume of the local heat conduction terms should 
equal the prescribed boundary heating rate. The sum over 
the volume of the local convection of energy or momentum 
terms should equal the boundary values. The sum of the 
cell densities divided by the total number of cells should 
equal the average density computed from the initial den- 
sity and the time integral of the mass density loss due to 
fluid expulsion. The change in global fluid pressure 
should equal the potential pressure decay when the non- 
isothermal cryogen is adiabatically restored to an iso- 
thermal state. The adiabatic potential pressure decay is 
defined algebraically as 

AP U* ■ («> < 24) 

£ v [lf >e > l 5P jS - ( -P e >ij] 

where <J > p is 0 at the average density and initial pressure, 
V is the cell volume, V c is the container volume, (EF) is 
the elastic factor of the system, Q DA is the adiabatic heat 
deficiency, (/)e) p ,9 is t ^ ie Product of the average 
density and the specific internal energy which is evaluated 
at the average density and initial global pressure, and 
( p e) ij is the product of the cell density and specific 
internal energy. The validity of equation (24) is pre- 
dicated upon the invariance of <J> over the pressure change 
( A P j max ) . 

Single-phase helium, parahydrogen, nitrogen, and oxygen in 
the density region from about the solid density through the 
low gas density exhibit a very weak or no dependency of <$> 
upon Pg at any given density. Thus, equation (24) is a 
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simple means of computing the adiabatic potential pressure 
decay. Where tank wall or heater tube, thermal capacitance 
is involved, and/or where boundary heating is prescribed, a 
pressure decay however abrupt may not be adiabatic. If the 
end state temperatures of the metal are known, the energy 
added to the system can be computed and for such a condi- 
tion equation (24) may be modified as 



Qdt 


UF) 


Q D T s - L V [<P e ><d, P/5 - (25) 

where Qb is the total boundary heating and Q M is the 
total energy from the metal for any short interval of 
time. Short here is taken to mean a time period in which 
the change of p is small enough so that the desired ac- 
curacy of equation (25) is maintained during the event. 


If the GEM is properly coded, all of the integral relations 
are satisfied exactly except for the energy convection 
terms which have a small nonconservative error. There is 
also a time dependent truncation error. In spite of these 
errors, the integral equation (24) and (25) are very nearly 
satisfied during time dependent calculations which is shown 
in an example in the results section. 

The 1/2 Box Model 


In order to simulate a physical problem with the GEM cer- 
tain boundary and initial conditions must be established. 
Figure (2a) shows a sketch of the cross section of a flight 
oxygen tank in the plane of the "g" vector which is taken 
to be perpendicular to the heater tube. Figure (2b) shows 
a sketch of the geometry used for the numerical model. The 
region shown has a depth dimension so that the walls of the 
model form a box. The left wall is assumed to be a plane 
of symmetry in order to reduce the total volume of region 
of interest by two. (This reduces the computer time by 
a factor of two) . For this reason, the model is called 
the 1/2 box model. Only one half the values of the con- 
tainer volume, heat input, expulsion rate, the heater tube 
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thermal capacitance (HTC) , and plumbing volume are used in 
the calculations. The boundary conditions for the 1/2 box 
model are 

1) freeslip left wall except where the heater element is 

2) no-slip right, top, bottom and heater surfaces 

3) uniform heat flux at the top, right and bottom wall 
for the heat leak simulation 

4) insulated left wall except where the heater surface is 

5) a flat plate heater surface midway up to the left 
wall, parallel to the left wall and extending the 
distance of the depth dimension 

The initial conditions for the velocity and temperature 
fields, the global pressure, and the initial density, are 
chosen to be as consistent as possible with the conditions 
in the flight hardware for each simulation. All but the 
first two conditions are usually known accurately. Special 
measures must be taken to estimate these two initial con- 
ditions and these are discussed in the results section. 

The Physical Description of the Apollo 
Single-Phase Oxygen System 

Figure (3) shows a schematic of an Apollo single-phase 
oxygen storage and supply system. The oxygen storage con- 
tainer is essentially a thin wall pressure vessel covered 
by a super insulation system. The container is vented and 
filled through pipes which are imbedded inside the super 
insulation for some distance before penetrating either 
the inner or outer container shells. On leaving the con- 
tainer, these pipes pass to the space craft surface. They 
are only active during the preflight oxygen loading pro- 
cedures and then are sealed off shortly before the flight. 

A third line is attached to the tank and supplies oxygen 
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to the fuel cell regulator and the surge tank. The func- 
tion of the surge tank is to separate the crew compartment 
from the low temperature cryogen and provide a large source 
of amb ient gas ready for use at any time. 

Inside the oxygen container are mounted two perforated 
tubes. One tube is a quantity gauge. The other tube is 
an extended surface for heat transfer. Several electrical 
resistance heater elements are helically wound around and 
welded to the heater tube. The control of the electrical 
power to the heating elements is accomplished by switches 
which are activated by fluid pressure. If the pressure 
falls below about 860 psia, the electrical power is engaged 
If the pressure rises above about 900 psia, the electrical 
power is disengaged. Prior to Apollo 14, two fans were 
installed in the heater tube, one in each end. When oper- 
ating, these fans draw oxygen through holes in the center 
of the heater tube assembly and ejected it radially at the 
ends, thereby stimulating circulating of the bulk of the 
stored cryogen. These fluid motions are potent enough to 
eliminate significant fluid temperature variations which 
are induced by heat transfer. As a result, the pressure 
history of the stored cryogen may be predicted by equa- 
tions (A10) and All). If the fans are not operated, the 
detail fluid motions in the container oxygen must be 
modelled in order to predict the pressure history. Success- 
ful modeling however, depends upon the prediction of the 
time variant "g" level, expulsion rate, and heat input. 

The method of determining each of these is now presented. 

Reduction of Apollo Flight Data 

It is the purpose of this discussion to provide some under- 
standing of the considerations which are made in deter- 
mining the "g" level, oxygen flow rate, and heat input 
histories as a function of the average stored oxygen 
density. 

Figure (4a) shows a sketch of the space craft in relation 
to the roll, pitch and yaw axis and the earth. For simpli- 
city the moon is not shown. Also shown in Figure (4b) is 
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a cross section of the space craft perpendicular to the 
roll axis and through the two Apollo 12 oxygen tanks. The 
heater tube in each tank is parallel to the roll axis. The 
"g" level and direction is defined if the inertial space 
acceleration vector can be evaluated. The method of rigid 
body dynamics is used to calculate the history of this 
vector. The angular rotation rates are taken directly from 
the Digital Auto Pilot (DAP) data. The angular accelera- 
tion is obtained by numerical differentiation of the 
rotation rates. When DAP data are not available, rotation 
and acceleration vectors are obtained by successive dif- 
ferentiation of the gyroscopic gimbal angles. The cal- 
culations showed that the accelerations induced by the 
solar pressure and the solar wind are orders of magnitude 
less than the accelerations induced by the thrust of noz- 
zles and rotations of the space craft about the three 
axes. (The Coriolis force is sometimes the same order of 
magnitude as the centripetal acceleration but for simpli- 
city this component is neglected for the 1/2 box model 
simulation. Variations of the acceleration magnitude with 
position in the stored fluid is also neglected for simpli- 
city in the 1/2 box model simulation.) 

There are five predominate types of flight modes which 
each have a typical "g" level and direction. They include 

1) engine burns 

2) attitude control maneuvers 

3) passive thermal control (PTC) 

4) attitude hold in translunar or transearth orbit 

5) attitude hold in lunar orbit 

The acceleration vector is parallel to the heater tube for 
type 1) operation. For the last three types of operation, 
the acceleration vector is predominantly perpendicular to 
the heater tube. Typical "g" levels in the last three are 
respectively 3 x 10“° "g" , 7 x 10” 8 "g", and 5 x 10“ 7 "g". 
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Abrupt increases in "g" level ("g" spikes) can accompany 
the use of the attitude control system. Then these "g" 
directions are not necessarily perpendicular or parallel 
to the heater tube. 

The certainty with which the " g" level is known is esti- 
mated to be better than a factor of 2 at 7 x 10”® "g" and 
better than a factor of 1.2 at 3 x 10”® "g". 

The energy input to the oxygen tank is composed of the 
heater power and insulation heat leak. The latter never 
exceeds about ten percent of the former. The estimation 

of the heat leak is no more accurate than the estimate of 
the oxygen flow rate but because it is a small fraction 
of the total energy input and because the heater power is 
known to about two percent, the overall energy input is 
believed to be known to better than three percent. 

The average cryogen density is obtained from the quantity 
gauge and time integrals of the expulsion rate, initial 
density, and container volume. The quantity is believed 
to be known to about two percent. 

The oxygen pressure is measured by a sensor whose output 
is digitized in four psi intervals at the rate of about 
one per second. The absolute value of pressure is known 
to about 20 psi at 900 psia. The differences in pressure 
are known to better than about four psi out of forty. A 
similar situation exists in the absolute pressure measure- 
ment of the surge tank. 

If the difference of the container and surge tank pressure 
depended only on the difference in these absolute pressure 
measurements, considerable error would be possible. How- 
ever, with normal crew compartment oxygen consumption 
rates and no use of surge tank gas for other purposes, 
the surge tank pressure must be the mean of the limit cycle 
values of the container. 

Thus a reference value of pressure is known and the surge 
tank pressure gauge can be scaled accordingly. The flow 
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through the restrictor depends upon the pressure differ- 
ential between the container and the surge tank. Since 
there is a flow meter between the crew compartment and 
the surge tank, a flow calibration of the restrictor is 
possible. The only time that the measurement of the flow 
to the surge tank requires the use of the restrictor cal- 
culation is when the surge tank pressure is different than 
the mean of the limit cycle values for the container. The 
estimate of the flow measurement accuracy is about 30 per- 
cent for the restrictor calculations. The mass flow rate 
to the fuel cells is known to about three percent which 
frequently is over half of the flow demand. The precision 
of the flow measurements then can be anything between 
three and thirty percent accurate depending upon the 
steadiness of the pressure in the surge tank and the crew 
compartment consumption rate. 

The only other way to predict the average flow to the surge 
tank is with equations (A10) and (All) . The heat input, 
cryogen density, and changes in cryogen pressure can be 
used to estimate the average flow rate to the surge tank. 
However this procedure only is adequate when the stored 
cryogen is very recently stirred by the fans so that it 
is spatially isothermal. 

It is hoped that the previous comments give some idea of 
the difficulty in obtaining accurate data which may be used 
to make predictions with the 1/2 box model. In spite of 
this situation, the predicted oxygen tank pressure history 
agrees very favorably with the flight data. Because 

of the data uncertainty in some cases, good agreement may 
just be fortuitous. Some reduction of the ambiguity can 
be made by parametrically varying various of the data 
which are most suspect (e.g., "g" level and flow rate). 

The trends of such analysis can be used to gain confidence 
in the results and thus reduce the uncertainty. 
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The Comparison of Predictions and Flight Data 


Two sets of data from the Apollo 12 flight have been re- 
duced and summarized. Predictions for these cases have 
been made with the 1/2 box model. The comparisons are 
presented after some introductory discussion, 

P s /P e is defined as the ratio of the time derivative of the 
global pressure in cryogen which has a nonuniform and uni- 
form time dependent temperature field respectively. For an 
ideal gas in which the compressibility and specific heats 
are constant, Ps/Pe is unity whether or not temperature 
variations exist. In nonisothermal cryogen, at a density 
of about twice the critical, P s /P e can approach a maximum 
value of about 16. As the density is reduced , P s /P g 

approaches a value of unity. (Note that if P s /P e is not 
unity, AP| max can acquire nonzero values.) The actual value 
that P g /P e acquires in finite difference calculations de- 
pends upon the resolution of the thermal boundary layer. 

With increasing acceleration level, natural convective mass 
rates increase near the heater surface and the heater sur- 
face temperature decreases. P s /Pe tends to unity for such 
heater temperature variation irrespective of the density of 
the cryogen. Coarse grid calculations of a thermal boun- 
dary layer with the 1/2 box model provide less boundary 
layer temperature variation and thus P s /P e tends to be less 
than if fine grid calculations are used. Accordingly, pre- 
diction of the pressure history is very sensitive to ther- 
mal boundary layer resolution where P s / p e potentially can 
deviate from unity which is the situation near and above 
the critical density. 

At cryogen densities below and above the critical density, 
the accuracy of the predicted pressure history is sensi- 
tive to thermal boundary layer resolution but for another 
reason. The heater tube has significant thermal capaci- 
tance. The coarse grid calculations always predict exces- 
sive heater surface temperatures and as a result of the 
heater thermal capacitance excessively long pressurization 
cycles are predicted. Convergence studies are employed to 
estimate the effects of coarse grid calculations. This 
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means that the same problem is repeated for a number of 
different grid sizes. 

For example, a 10 x 10, 20 x 10, 30 x 10, 40 x 10, 60 x 10, 
80 x 10, 20 x 20 (the first and second integers are the 
number of computational cells in the x- and y- directions 
respectively), etc., grids may be used to predict a pro- 
blem. By appropriately plotting the data, the convergence 
characteristics of some parameter of interest can be 
determined. If this procedure is used in a careful way 
beginning with the coarse grids first, estimates can be 
made as to how fine the grid must be to get the desired 
accuracy. This will be illustrated later. 

The potential pressure decay (AP|max) ma y k e defined as 
the amount the global pressure should decay when a 
cryogen is adiabatically restored to an isothermal state 
during which the fluid expulsion rate is zero. Trunca- 
tion and conservation of energy errors can be delineated 
if the time integral of Pg is compared to the difference 
in the state of AP| max for the boundary conditions of no 
boundary heating, no fluid expulsion, and no heater thermal 
capacitance. To illustrate such a comparison, two problems 
have been chosen and they differ only in the grid sizes 
which are a 6 x 6 and 12 x 12 and in the time step used 
after the "g" spike is initiated. The box has the dimen- 
sions of 2 x 2 x 1.187 feet for the height, width, and 
depth respectively. The initial density is about twice the 
critical density. Figure (5) shows the history of Pg and 
A P|max °f the 12 x 12 grid problem in which the fluid 
expulsion, wall heating, heater input, and acceleration 
level are nonzero for a period of time. The wall heating 
and heater input are chosen in order to generate an oscil- 
lating pressure history. The heater tube thermal capa- 
citance is zero during the entire duration of the problem. 
The "g" level is chosen at 2 x 10” 8 "g" so that the buoy- 
antly driven mass rates are of the same order as the expul- 
sion rate. At some later time, the "g" level is increased 
to 2 x 10~ 8 "g" in a time step and at the same time the 
expulsion rate and the boundary heating are set to zero. 

The "g" level increase causes a substantial increase in 
the buoyantly driven fluid velocity. As a result, the 



high energy fluid which is concentrated in a glob near the 
heater surface is rapidly stretched which substantially 
enlarges it's surface area. This in turn enhances the 
heat conduction and results in a rapid decrease in fluid 
temperature differences. The sharp reduction in nonuniform 
temperature field causes the global pressure to sharply 
decay. The pressure decay rate is an order of magnitude 
sharper than that which would have occurred had the "g" 
level been unchanged during the course of the problem. 

The time step is 30 seconds prior to the "g" spike and it 
is reduce d by 
n , / (AX, 

At = 0.4 y 

when the "g" spike occurs which yields a A t = 0.6977 se- 
conds and At = 0.4933 seconds for the respective problems. 
The change in P g of 41.0 and 49.3 corresponds to a change 
in A P | max of 42.5 and 49.6 respectively for the 6x6 
and 12 x 12 grid problems. This yields respective differ- 
ences in the actual and ideal decay changes of 2.8 and 
0.3 psi for the respective problems. The reduction in 
the error is entirely due to the time step reduction. Fur- 
ther reduction in the time step after the "g" spike will 
not significantly improve the result. This is due to the 
fact that there is a small time dependent nonconservative 
error in the energy convection terms. This error can be 
virtually eliminated in the manner outlined in reference 
( 2 ) . 

Finally, the GEM has one other time related truncation 
error which for the problems considered is trivial. The 
algorithm exhibits a slight damping of kinetic energy and 
this error may be eliminated with a nontrivial increase 
in computer time. This consideration will be left to a 
future paper. 

As a rule, the truncation errors in time are well control- 
led as the grid is refined. The significant errors for 
fine grid calculations reduce to the definition of the 
initial and boundary conditions. As has been mentioned 
previously the mass expulsion rate is probably the most 


75 



inexactly known quantity® There are some other significant 
ones which are now discussed. 

The elasticity factor (EF) for which a derivation is provi- 
ded in the Appendix can be evaluated precisely except for 
term ( p V n /P c V P ) which accounts for the exter- 
nal volumli ef¥ect. Since there is no temperature sensor 
anywhere on the external plumbing, estimates of p and n 
(which is the reciprocal of the polytropic coeff?cient) 
must be made. (It is believed that p p should be about 
5 lb^ /ft and n about 1.0. In this case, the magnitude 
of the effect of the external plumbing upon the elasticity 
factor is nearly negligible.) The elasticity factor is 
about one half at the high densities and rapidly approaches 
unity as the density is reduced to the critical value. It 
is unity for densities less than the critical value as well. 
Another difficulty in using the 1/2 box model arises in 
choosing the initial conditions for a given problem. As a 
rule, a uniform initial temperature field with a zero velo- 
city field are the chosen initial conditions. The problems 
selected for analysis always have a fan cycle prior. Be- 
cause the pressure response of the system appears to satis- 
fy equations (A10) and (All) of the Appendix, the isothermal 
assumptions appear to be warranted. But the zero velocity 
field prescription is certainly wrong except after the fluid 
motions from the fan cycle have decayed below those gener- 
ated by the buoyant forces. Estimates of this decay time 
are made with fine grid calculations in which the initial 
field values of velocity are large. It has been found 
that less than about two hours are required to meet this 
condition at all fluid densities. 

The treatment of the heater tube as a vertical flat plate 
in the 1/2 box model is believed to be wrong in one res- 
pect. Since the heater tube has only a few square inches 
of hole area through which high and low density fluid may 
flow, at low enough "g" levels some heated fluid can be 
trapped inside. This has the net effect of reducing the 
effective surface area for heat transfer. The "g" level 
at which the effective surface area is significantly re- 
duced is at about 10“^ "g" or less. 
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The extremes of the surface area possible are 

1) the total of the outside and the inside of the heater 
tube (about 1.1 ft 2 ) 

2) the exterior of the portion of the tube which is 
wrapped with the heating elements (about 0.475 ft 2 ) 

o 

For the calculations presented, the value of 0.95 ft is 
used and is twice the value of 2) . This area certainly is 
a bit too large after some time period because of the 
restraint which the tube holes impose on the buoyantly 
driven flow and because the heater tube is not exactly 
isothermal along its length. (Evaluation of these two 
effects is currently under study.) 

With these various reservations in mind, the two problems 
are now considered. 

Two problems from Apollo 12 flight data were chosen for 
analysis because initial and boundary conditions could 
be reasonably assessed and because these problems exhibit 
the peculiar effects that the nonisothermal cryogen has 
upon the pressure history. 

The 1/2 box model is assumed to be one foot wide, two feet 
high, and the depth is computed so that the volume is that 
of 1/2 of the oxygen tank which is 4.75 cubic feet. 

The first flight problem of interest involves the predic- 
tion of the time required to raise the oxygen pressure 
about 33.5 psi in twelve minutes and forty seconds. The 
fan cycle occurred about eight hours prior to the time of 
interest. The pressure decayed steadily and the expulsion 
rate history could be found from equations (A10) and (All) . 
These values agreed well with those predicted from crew 
compartment and fuel cell flows. For this reason and since 
only heat leak was involved during the decay cycle, the 
system was believed to be essentially isothermal when the 
heater was activated, which is the beginning of the pro- 
blem. The eight hour period was more than enough time for 
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the fluid velocity field to decay adequately so that the 
zero field velocity prescription could be used at the be- 
ginning of the prediction with the 1/2 box model. The 
elasticity factor is about 0.97 so that the question of 
how to evaluate the external volume effect is trivial. 
Figure (6) shows the results of the 1/2 box model predic- 
tions as a function of the grid size with and without the 
heat tube capacitance of 0.07 BTU per degree Rankine. It 
may be noted that the 60 x 10 and the 80 x 10 grid results 
are nearly identical. Comparison of the velocity profiles 
show that the boundary layer is well resolved by the 80 x 
10 grid. 

It is curious that the resolution is only weakly dependent 
upon the number of cells in the y- direction. Figure (7) 
shows the asymtotic estimate of the 1/2 box model predic- 
tions. Figure (8) shows these asymtotic values replotted. 
The agreement of the prediction and the flight data is 
within about two percent which is better than would be 
expected in view of the probable errors of the various 
measurements. There are some interesting things to 
note about the predictions. 

The cryogenic oxygen density is slightly less than the 
critical density. The average rate change of pressure is 
about 15 percent larger than that which would occur for 
isothermal cryogen. (P s /Pe is about 1.15 average value. 

As a result, the adiabatic potential pressure decay ac- 
quires a value of about 4 psi during the pressurization 
stroke.) Also, the heater tube thermal capacitance is 
shown to be a significant factor in the length of the 
pressurization stroke. Thus without accounting for the 
transient heater tube and thermal boundary layer inter- 
action, the flight data could not well be predicted. 

The 1/2 box model is used to simulate this problem again 
but it differs from that shown in Figure (6) only in the 
"g" level which is set at 4 x 10“^ "g". This lower "g" 
level results in lower buoyantly driven mass rates by the 
heater surface which results in a larger asymptotic 
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heater surface temperature. The asymptotic temperature 
differences between the heater surface and the bulk fluid 
is 230° R and 340° R for the problems at 2 x 10“ 6 "g" and 
4 x 10 “^ "g" respectively. Had the P s /P e ratio remained 
the same for each of these problems, the latter problem 
would have had a pressurization stroke about a minute 
longer than that found for the former problem. This time 
increase is due entirely to the increased energy require- 
ment of the heater tube. However, since the average value 
of P s /P e for the latter problem is 1.23, the time increase 
caused by the heater tube capacitance effect is almost 
exactly offset and the resultant stroke length is only a 
few seconds longer. Since the reduction of the heater 
surface area to about a value of 0.7 ft^ in the former 
problem would result in about the same asymptotic heater 
surface temperature, the pressurization stroke length would 
be almost exactly the same as that found in the latter pro- 
blem. But since the latter and former problems yield al- 
most identical stroke lengths, this problem is not suited 
to determining with precision the effective heater tube 
area. After examining the various errors and their effect 
on the results of the simulation, it is concluded that the 
heat tube thermal capacitance and the differential pres- 
sure measurement during the stroke are the primary para- 
meters upon which the accuracy of the simulation depends. 
Since the agreement is excellent between the flight data 
and the prediction, it is concluded that whatever errors 
there are in these two parameters, they must be compen- 
satory. 

The second flight problem of interest examined in detail 
involves a high density problem in which the initial den- 
sity is about twice the critical density. The elasticity 
factor is between 0.58 and 0.35 depending upon how the 
external volume effect is accounted for. The latter value 
assumes that the entire external plumbing volume is filled 
with cryogen which is at the density of the average of that 
of the stored fluid. From thermal studies of this plumbing 
system, such a condition is unlikely and it is probably 
the case that the entire volume contains near ambient 
oxygen. Since the process in the plumbing during the 
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pressure oscillations is probably isothermal, n is probably 
unity. With these assumptions the elastic factor is about 
0.55 which means that the contribution from the external 
volume is only about five percent of the total. 

The flow rate is time dependent and is modeled with the 
expression 

m ■ (0.755 + A + B) 


where 


A = 0.1 lb^/hr when the heater is on 


A = 0 lb m /hr 
B = 0.17 (P g 
B = 0 if (P g 


when the heater is off 
- P« 


. 0.488 , 

surge' J -- L vr g 


P surge^ ^ ® 


p surge) > 0 


where B is the flow through the orifice which is upstream 
of the surge tank. B is zero if the pressure difference 
(Pg - Psurge) is negative because a check valve in the line 
prevents reverse flow. P SU rge is the surge tank pressure. 
The "g" level is time dependent and the values used in the 
simulation are shown in Figures (9, 10, and 11). Actually 
there are many more "g" spikes than indicated but these 
were deemed insignificant compared to those shown. Figures 
(9, 10, and 11) also show the end points of the pressure 
strokes for the flight data. The fan cycle occurred at 
about a ground elapse time (GET) of 4 hours and 27 minutes 
so the residual velocities should have been gone at the 
most by the time GET = 6:27. However examination of the 
length of the strokes show that within about 40 minutes the 
predominate effects of the fan cycle had dissipated. Thus 
the stroke length is then relatively fixed and is modified 
predominantly by the "g" spikes. Note how sharp the pres- 
sure change is with the "g" spikes but that no substantial 
change in fluid pressure occurs until the "g" level in- 
creased to 10~3 "g" . a dramatic decay of fluid pressure 
then ensues and an enlarged view of this event is shown in 
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Figure (11) . The fans are turned on before the course of 
the "g" spike induced pressure decay stroke is run. It 
appears however from the slope of the decay curve prior to 
the fan cycle, that virtually all of the potential for 
pressure decay would have been absorbed. Thus properly 
interspersed "g" spikes of sufficient magnitude could be 
used to provide the same function as the fans, that of 
stirring the stored cryogen so that it's temperature field 
is for practical purposes isothermal. 

Attempts to simulate the pressure history of this problem 
have begun but are not yet complete. As a consequence, 
only preliminary results from 1/2 box model studies are 
presented. Between GET = 4:29 and 4:41, the decay cycle 
as noted before, involves a cryogen whose temperature field 
is isothermal, but time dependent. The end points of this 
pressure cycle have been predicted to better than twenty 
seconds in eleven minutes with the 1/2 hox model. The same 
result is obtained with the equations (A10 and All) . So 
that the remainder of the four hour simulation could be put 
on one plot, this decay stroke is omitted from any display. 
Figure (12) shows a pressure trace for a 60 x 10 grid in 
which the previously given flow rate equation was not em- 
ployed. (The pressure trace for the included flow equa- 
tion could not be prepared in time for this publication. 

With some exceptions the average flows are about the same 
and for discussion purposes it will serve to illustrate 
the important features of the simulation.) A convergence 
study shows that the 60 x 10 grid almost adequately resolves 
the thermal boundary layer during the periods between GET = 
4:41 and 5:53 and between 6:40 and 8:33. During, and some 
time after the "g" spikes, the thermal boundary layer is 
not well resolved. It is estimated that about at least a 
150 x 40 grid would be required during this period. (With 
existing computers, this degree of resolution is impracti- 
cal. However, a variable mesh grid system is being de- 
veloped which hopefully will permit resolution which is 
more nearly adequate.) 

It should be noted that the length of the pressurization 
stroke in Figure (12) is short sooner than shown in Figure 
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(6) . This is believed to be due to the difference in the 
initial conditions which has been already noted. At GET 
5:41, the stroke lengths are about the same and this may 
be noted by examining Figure (13) which shows the plot of 
the pressurization stroke length versus the number of grid 
points. (The flight data is also shown.) . The asymptotic 
value is believed to be almost attained with the 80 x 10 
grid. (The variable mesh program may provide the tool 
necessary for the practical verification of this belief.) 

If so, the agreement between the simulation and the flight 
data is good. 

The stroke length of the pressurization cycle shown in 
Figure (13) is about three times less than that predicted 
for a cryogen whose temperatire field is isothermal, but 
time dependent. If the heater tube thermal capacitance is 
zero in the 1/2 box model, the stroke length is about 
fifteen times less than that predicted for a cryogen whose 
temperature field is isothermal, but time dependent. The 
elasticity factor and expulsion mass rate inaccuracies 
can, at most, only affect the pressurization stroke length 
by about a factor of two. Thus accurate simulation of the 
transient thermal boundary layer and heater tube tempera- 
ture are the most important considerations in predicting 
the pressurization stroke length. 

With regard to the simulation of the "g" spikes, it may be 
noted that there are several qualitative and quantitative 
resemblences with the flight data. This is true of several 
of the small "g" spikes and the large "g" spike event which 
occurs at the end of the simulation. This latter event has 
been enlarged and displayed in Figure (14) . The magnitude 
of the pressure decay is substantially attenuated by the 
heat which is derived from the rapid cooling of the heater 
tube. (Convergence studies are currently being performed 
for this "g" spike event in order to attempt to gain quan- 
titative data on the grid sizes required for adequate 
resolution.) The obvious disagreements between the flight 
data and the simulation are believed to be due to 

1) since simulated strokes are not in phase with flight 
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"g" data, the "g" spikes occur at different relative 
positions in the stroke than is shown for the flight 
data 

2) the geometry of the heater tube probably plays some 
part 

3) the grid is certainly too coarse 

Further studies of this problem will hopefully provide the 
information which is necessary to delineate the importance 
of each of these considerations. 

In spite of the preliminary nature of the data thus far 
acquired for this flight problem, the results to date are 
very encouraging, enough so that further work is warranted. 

Conclusions 


The methodology of the General Elliptic Method for the simu- 
lation of time dependent single-phase cryogenic flow has 
been presented. A brief description of the Apollo 12 sin- 
gle-phase oxygen storage system has been presented along 
with some discussion of the elements which are considered 
for the reduction of flight data. The 1/2 box model in 
which the boundary conditions used to simulate the cryogenic 
oxygen flow processes of the Apollo 12 oxygen tank has been 
presented. Mention has been made of some of the difficul- 
ties which have arisen in obtaining accurate quanitative 
data from the 1/2 box model and the flight data. In spite 
of these difficulties, the comparisons of the predictions 
and the flight data are shown to be good. 

It is recommended that further application of the 1/2 box 
model be made to predict flight system response as the 
flight data becomes available. Assessment of the deficien- 
cies and capabilities of the 1/2 box model would then be 
possible for a wide range of flight conditions. If the 
overall results are reasonably favorable, prediction of the 
oxygen tank pressure history for future missions could then 
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be undertaken with a confidence well beyond that pre- 
viously available. 
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Figure 1. TYPICAL CROSS SECTION OF A COMPUTATIONAL CELL 
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Figure 3. APOLLO 12 OXYGEN TANK AND SUPPLY SYSTEM SCHEMATIC 
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ACCELERATION DUE TO ROTATION: 

A rot -(&xR) + (WxWx R) 

ACCELERATION DUE TO GRAVITY GRADIENT: 

A * * 

-2M e GR e «R * e 

a gg * a total = a rot + a gg 
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SOLAR PRESSURE AND SOLAR WIND ARE NEGLECTED 

Figure 4. O 2 TANK ACCELERATIONS 
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FLUID PRESSURE (PSIA) HEATER TUBE TEMPERATURE ( R) 



NOTE: 

QUANTITY 35% 

HEATER POWER 
114 WATTS 

HEAT LEAK 
20 BTU/HR 

FLOW RATE 

1 LB/HR AT 860 
PSIA AND IS 
INCREASED BY 
0.37 (P g - 860)/ 25 

ACCELERATION 
2x10 6 

ELASTIC FACTOR 
0.975 

HTC - 0.07 - HEATER 
THERMAL 
CAPACITANCE 
RISE TIME OF 

p ref + ap f“ 

12 MIN 40 SEC •) 


REF 


• 859.8 
33.5 


20 


TIME (MIN) 


Figure 6. APOLLO 12 OXYGEN TANK No. 1 
PASSIVE THERMAL CONTROL MODE 
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MAX POTENTIAL PRESSURE DECAY (PSD 




TIME (MIN) HEATER SURFACE TEMPERATURE (°R) 



20 40 6 0 80 100 120 


NUMBER OF x- DIRECTION CELLS 

Figure 7. APOLLO 12 OXYGEN TANK No. 1 
PASSIVE THERMAL CONTROL MODE 
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AVERAGE TUBE TEMPERATURE, T (®R) 


23 PROBABLE ERROR BETWEEN 
PREDICTION AND EXPERIMENT 


T (OR) 


r AP (PSIA) 


NOTE: 

QUANTITY 35% 

HEATER POWER 114 WATTS 
HEAT LEAK 20 BTU/HR 
FLOW RATE 1 LB/HR AT 860 
PSIA AND IS INCREASED 
BY 0.37 (P g - 860/25) 

ACCELERATION 2 x 10" 6 
CONTROL 859.8 AND 900 PSIA 
ELASTIC FACTOR 0.975 
o FLIGHT DATA 
A P - 33.5 PSIA 
TIME -12 MIN 40 SEC 
P s /P e a 1.14 


TIME (MIN) 


Figure 8. APOLLO 12 OXYGEN TANK No. 1 
PASSIVE THERMAL CONTROL MODE 
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OXYGEN TANK ACCELERATION OXYGEN TANK PRESSURE 

LOG 10 V <PSIA) 




Figure 9. APOLLO 12 FLIGHT DATA - ATTITUDE HOLD 
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Figure 10. APOLLO 12 FLIGHT DATA - ATTITUDE HOLD 



OXYGEN TANK ACCELERATION OXYGEN TANK PRESSURE (PSIA) 
LOG™ g 






QUANTITY 95 PERCENT HEATER POWER 114 WATTS 
ACCELERATION 7X10 -8G HEATER AREA 0.95 FT. SO. 
GRID SIZE 60 X 10 THERMAL MASS 0.07 BTU/DEG R 



Figure 12. APOLLO 12 SIMULATION 
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TANK PRESSURE (PSIA) 



GROUND ELAPSED TIME, GET (HOURS: MIN) 


Figure 13. APOLLO 12 "g" SPIKE SIMULATION 
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NOTE: 

QUANTITY 95% 

HEATER POWER 114 WATTS 
HEAT LEAK 14.5 BTU/HR 
FLOW SHOWN IN FIGURE 13 

ACCELERATION 7 x 10“ 8 
CONTROL 860 TO 900 PSIA 
ELASTIC FACTOR 0.397 
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Figure 14. APOLLO PRESSURE STROKE RATE VS GRID SIZE 
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Appendix A 


The Derivation of the Elastic Factor (EF) for the 
Spherical Container and the External Plumbing 



System Schematic 


Assumptions: 


1) The volume (V c ) is enclosed by a thin walled elas- 
tic spherical container with a response of 
/ J_\ 4YI - 3r(l“ _gj 

where 


Zbt 


r = container radius 
cr = Poisson's ratio 
b = container thickness 
E = Young's modulus 

2) The plumbing volumes (V p ) are rigid (dV/dP)| p = 0 and 
contain a compressible gas with a response P = 

Z p RT in which the rate change of compressibility 
(Z) with respect to the pressure (P) is zero 
(dZ/dP = 0) . 

3) The time (t) derivative of pressure is a constant 
for all the fluid in the system (dP/dt)| c = 

(dP/dt) / p . 
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4) The enthalpy and temperature properties are only 
functions of the pressure and the density for the 
stored fluid. 

Now define 


e 

_aWl 

* 

« »P-1 , 

*(p«) Ip 

h s 

e ♦ £ 

dMc 

dt 

s - iric * 

drt p 

at 

• • 

s m c - rop j 

dP _ 

at * 

• 

P 

dP 

It 

• 

« P 1 

/°c 

* (.-v) c ’ 

<°p a 

t^-) p 


The result of time differentiation pf^and/Op after in- 
voking assumption 2) and some appropriation definitions 
are 


• = dM I Mr dVo _ P C dv c 

v c at | c " at “ v c v c at 


(Al) 



(A2) 


The energy balance for V c is 


d_ 

at 


(P e v) c 



m c h c - P 


av c 

dt 


(A3) 
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The subscript c is now dropped from (A3) until some inter- 
mediate manipulations are complete. Since 


/oe = /o e (p, p) 

- + 1 -^ h - e) + l (A4) 


Now by substitution eliminate d (^ej/dt from (A3) and 
(A4) and combine the result with (Al) and an appropriate 
definition to yield 





dv 

dt 


C/O0) 


a-me 


(A5) 


Now by the chain rule 


dV__dV dP 
dt " dP dt 


and combined with (A5) with the subscript c restored is 


P c - (-f) c (ic-^c 9 C>UF) 

where 

(tF) * 

Invoking assumption (1) , (A7) may be written as 


(A6) 


(A7) 
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(A8) 


UF) 


\_ 

\ + (,p<t>0) c 


3r (i -cr) 

2bE 


The state equation P = ZR|OT with the appropriate defini- 
tion may be written as 


p » (ze «l) f 

and may be differentiated with respect to Vp and after 
some alegraic manipulation is 



and after assumptions (3) and (4) and the definition of 
dpp/dt is employed, simplifies to 


• (m c - fifip) Pp 

P P " (<oV) p 

(A6) , (A8) , and (A9) may be combined to yield 

P c - <■$),, t «c> <•«=> 


(A9) 


(A10) 


where 
(tF) = 


l + 


(p4>e) c 


\ 

3r (i-o*) 
Z.bE 


where n = 1 


(pV) p (")‘ 
Cp>»)c p > ■ 


(All) 


It can additionally be shown that for a polytropic heating 
1 ^ n ^ *- v /C where n is unity for the isothermal case and 
C v /Cp for thfe isentropic limit. 
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